Note
This example notebook demonstrates NSphere visualization techniques. The original notebook can be found in the examples directory of the NSphere project.
Copyright 2025 Marc Kamionkowski and Kris Sigurdson. Licensed under Apache-2.0 (http://www.apache.org/licenses/LICENSE-2.0).
This notebook provides some simple python scripts to make plots and movies from the output of NSphere-SIDM.
The scripts use data produced by nsphere --sidm --ntimesteps 12500000 --dtwrite 25000 --nout 250 --tfinal 1100 --sort 3 (the parameters used for the plots in the paper). This is then run with 100000 particles, the default halo profile, an NFW profile with scale radius \(r_s = 1.18\) kpc, mass \(= 1.15 \times 10^{-9}\) M\(_\odot\), and a large-radius cutoff function \([1+(r/r_{\text{vir}})^{10}]^{-1}\) with \(r_{\text{vir}} = 19\) (concentration parameter). The SIDM
cross section (per unit mass) is isotropic and velocity-independent, and taken to be 50 cm\(^2\)/g. All of these can be changed using command line options, run nsphere --help for more information.
The jupyter notebook is set up to be run from the NSphere-SIDM directory. It was run with Python 3.12.4
SIDM Core Collapse Analysis, Visualization and Plots
[1]:
# Import necessary libraries
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from scipy.ndimage import gaussian_filter # 2D smoothing filter
from IPython.display import HTML
from matplotlib.colors import PowerNorm
import os
# Define the path to the data files and the time steps
# Adjust the path to your data files as necessary
path_template = "../data/Rank_Mass_Rad_VRad_unsorted_t00{timestep:03d}_100000_12500001_1100.dat"
time_steps = range(0, 501, 2) # Time steps (251 snapshots from 0 to 500)
# Convert time steps to actual time in Gyr
time_in_gyr = [19.605 * t / 500 for t in time_steps]
# Define bin edges for radius (R) and radial velocity (Vrad)
# Fine-tuned for optimal resolution while staying under 20MB limit
bin_edges_R = np.arange(0, 10, 10/154) # 154 bins
bin_edges_V = np.arange(-50, 50, 100/154) # 154 bins
print(f"Using {len(bin_edges_R)-1} x {len(bin_edges_V)-1} bins")
# Define the record dtype (as in your original code)
record_dtype = np.dtype([
('rank', np.int32),
('mass', np.float32),
('R', np.float32),
('Vrad', np.float32),
('PsiA', np.float32),
('E', np.float32),
('L', np.float32)
])
# Storage for 2D histograms
all_histograms = []
# Loop over time steps and process each file
for time_step in time_steps:
file_path = path_template.format(timestep=time_step)
try:
# Read the binary file
alldata = np.fromfile(file_path, dtype=record_dtype)
# Extract radius and radial velocity data
data_R = alldata['R']
data_Vrad = alldata['Vrad']/1.023e-3 # Convert to km/s
if np.isnan(data_R).any() or np.isnan(data_Vrad).any():
print(f"Warning: NaN values found in file at time step {time_step}")
if not (np.isfinite(data_R).all() and np.isfinite(data_Vrad).all()):
print(f"Warning: Non-finite values found in file at time step {time_step}")
# Compute 2D histogram: counts for each (R, Vrad) bin
hist, _, _ = np.histogram2d(data_R, data_Vrad, bins=[bin_edges_R, bin_edges_V])
# Optionally, apply 2D Gaussian smoothing
# Adjust sigma=(sigma_R, sigma_V) as appropriate
smoothed_hist = gaussian_filter(hist, sigma=(0.1, 0.1))
# Optionally, you could apply a normalization/scaling here if desired
all_histograms.append(smoothed_hist)
except Exception as e:
print(f"Error reading file {file_path}: {e}")
# Append a zero histogram if file is missing or cannot be read
all_histograms.append(np.zeros((len(bin_edges_R)-1, len(bin_edges_V)-1)))
# Convert list to numpy array for easier handling
all_histograms = np.array(all_histograms)
print("Histogram data shape:", all_histograms.shape) # (time_steps, n_R_bins, n_V_bins)
# Determine common color scale for the animation
vmin = 0
vmax = np.max(all_histograms) # Use data maximum across all frames
# Create the animation with power-law color scaling
fig, ax = plt.subplots(figsize=(8, 6))
# Note: We transpose the histogram so that the horizontal axis corresponds to R and vertical to Vrad.
# Use PowerNorm with gamma=1/4 for x^(1/4) scaling
# Use 'hot' colormap for better distinction of high-density regions (black->red->yellow->white)
im = ax.imshow(all_histograms[0].T, origin='lower',
extent=[bin_edges_R[0], bin_edges_R[-1], bin_edges_V[0], bin_edges_V[-1]],
aspect='auto', cmap='viridis', norm=PowerNorm(gamma=3/4, vmin=vmin, vmax=vmax))
ax.set_xlabel(r"$r\;\;$[$\mathrm{kpc}$]")
ax.set_ylabel(r"$v_r\;\;$[$\mathrm{km}/\mathrm{s}$]")
ax.set_title(r"$t = %.2f\;\;$[$\mathrm{Gyr}$]" % time_in_gyr[0])
fig.colorbar(im, ax=ax, label=r"Counts per bin")
def update(frame):
im.set_data(all_histograms[frame].T) # update image data
ax.set_title(r"$t = %.2f\;\;$[$\mathrm{Gyr}$]" % time_in_gyr[frame])
return im,
ani = FuncAnimation(fig, update, frames=len(time_steps), repeat=True, interval=200, blit=True)
# Uncomment to show initial static image
# plt.show()
# Uncomment this line to save the animation as a GIF
# ani.save("results/phasespace.gif", writer="pillow", fps=10)
# Close the figure to prevent static display
plt.close(fig)
# For displaying in a Jupyter Notebook:
HTML(ani.to_jshtml())
Using 154 x 153 bins
Histogram data shape: (251, 154, 153)
[1]: